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Abstract 

Drug-induced liver injury (DILI) is a major public health issue and of serious concern for 
the pharmaceutical industry. Early detection of signs of a drug's potential for DILI is vital 
for pharmaceutical companies' evaluation of new drugs. A combination of extreme values of 
liver specific variables indicate potential DILI (Hy's Law). We estimate the probability of 



severe DILI using the Heffernan and Tawn ( 2004[ ) conditional dependence model which arises 



naturally in applications where a multidimensional random variable is extreme in at least 
one component. We extend the current model by including the assumption of stochastically 
ordered survival curves for different doses in a Phase 3 study. 

Keywords: stochastic ordering; multivariate extremes; conditional dependence; drug toxicity; 
liver injury 

1 Introduction 

Drug-induced liver injury (DILI) is a major public health and industrial issue that has concerned 



clinicians for the past 50 years. The FDA (2008) reports that many drugs for a diverse range 



of diseases were either removed from the market or rejected at the pre-marketing stage because 
of severe DILI (e.g., iproniazid, ticrynafen, benoxaprofen, bromfenac, troglitazone, nefazodone, 
etc.). Therefore, signals of a drug's potential for DILI and early detection can help to improve 
the evaluation of drugs and aid pharmaceutical companies in their decision making. However, 



in most clinical trials of hepatotoxic drugs, evidence of hepatotoxicity is very rare and although 
the pattern of injury can vary, there are no pathogenomic findings that make diagnosis of DILI 
certain even upon liver biopsy. Therefore, the majority of drugs that have been withdrawn fall 
into the post-marketing category and it is believed that the frequency of severe DILI does not 
exceed 0.01%. 



Although the mechanism that causes DILI is not fully understood yet, the procedure at which 
its clinical assessment is performed stems from Zimmerman's observation that hepatocellular 
injury sufficient to impair bilirubin excretion is a revealing indicator of DILI (Zimmerman 1978, 
1999), also informally known as Hy's Law. In other words, a finding of alanine aminotransferase 
(ALT) elevation, usually substantial, seen concurrently with bilirubin (TBL) greater than twice 
the upper limit of the normal (ULN), identifies a drug likely to cause severe DILI (fatal or 
requiring transplant) at a rate roughly 1/10 the rate of Hy's Law cases. Moreover, these eleva- 
tions should not be attributed to any other cause of injury; such as other drugs and alkaline 
phosphatase (ALP) should not be greatly elevated so as to explain TBL's elevation. 



Southworth and Heffernan (20126) identified the assessment of DILI as an extreme value problem 



and using the Heffernan and Tawn (2004) modelling approach, analysed liver-related labora- 



tory data. The use of the Heffernan and Tawn (2004) model in this context is supported by the 
flexibility of the model to allow a broad class of dependence structures and the possibility to 
describe the probabilistic behaviour of a vector random variable which is extreme in at least one 
margin. Despite its strong modelling potential, complications in terms of parameter identifia- 
bility problems and invalid inferences are experienced with the original modelling procedure of 



Heffernan and Tawn (2004). Keef et al. (2012) provided missing constraints for the parameter 



space of the Heffernan and Tawn (2004) model that are aimed to overcome these complications. 



The data we consider in this study relates to observed liver-related variables from a sample of 
606 patients who were issued a drug that has been linked to liver injury in a Phase 3 clinical 



trial and can be found in Southworth and Heffernan (2012c); see also Southworth and Heffernan 



(2012&). The patients were categorised into 4 different dose levels in a randomised, parallel group. 
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double blind Phase 3 clinical study. Our objective in this paper is to extend Keef et al. (2012) 
constraints to bound any two conditional distribution functions and hence incorporate scientific 
knowledge and efficiency into the modelling procedure of DILI by imposing the assumption of 
stochastically ordered survival curves in the joint tail region of the ALT and TBL variables 
between different dose levels of a drug. Therefore, estimation of ALT and TBL under the dose 
ordering assumption is beneficial as it sharpens inference and removes variability that arises in 
small sample sizes in the joint tail region of ALT and TBL that identifies DILI. Our motivation 
to impose additional constraints stems from the fact that the probability of DILI is logically 
ordered between different dose levels when the drug is liver toxic and this feature is unlikely 
to be evident from the data when considering clinical trials with small sample sizes as this 
particular one. 



Table 1: Kendall's r estimates between ALT and TBL 





Baseline 


Post-baseline 


Residual 


Dose A 


0.13 


0.25 


0.12 


Dose B 


0.09 


0.12 


-0.01 


Dose C 


0.11 


0.10 


-0.02 


Dose D 


0.15 


0.18 


0.05 



For example, the two central columns of Table [T] show the estimated Kendall's rank correlation 
of ALT and TBL measured at the baseline and post-baseline periods. The last column shows 
the Kendall's r of the residuals of ALT and TBL that were obtained from the median regression 
of the log-post-baseline on the log-baseline variable, see also Section 5.1 Table [T] conveys this 
lack of ordering since the estimated dependence of ALT and TBL at dose A appears to be sig- 
nificantly different from zero and even higher than any other dose in post-baseline and residual 
scale. On the other hand, the published literature reports jaundice, hepatitis and similar symp- 



toms in approximately 1 out of 500 patients taking the dose D of this drug ( Southworth and 



Heffernan, 20126). We view such high dependencies in low doses as a by-product of sampling 



variability. 



The proposed methodology and data analysis of the paper are based on an asymptotically mo- 
tivated model of multivariate extreme value threshold model which is fitted to a fraction of the 
data. An alternative approach would be to model the joint distribution between the variables 
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using all the data through the use of empirically selected marginal and copula models (Joe 



1997, Nelsen, 2006). The former approach should exhibit less bias but larger variability than 



the latter, thus as the sample size increases the extreme value approach is likely to become 
more the efficient. The sample size in our study is probably about at the boundary where the 
extreme value methods have an advantage. Furthermore, even if the copula approach were to 
be adopted, the strategies developed here would still be relevant as the stochastic ordering issue 
would need addressing. 



The paper is organised as follows. Section [2] describes the conditional dependence model of 



Heffernan and Tawn (2004) and the constraints of Keef et al. (2012). The additional constraints 



based on the assumption of stochastically ordered conditional distributions are presented in 
Section [sj The effect of the constraints of Keef et al. ( 2012[ ) and this paper are assessed with 
a simulation study in Section |4j We illustrate in Section [5] the application of the additional 
constraints by analysing the multivariate extremes of ALT and TBL of the DILI data. 



2 Methodology 



2.1 Marginal transformation 

Here and throughout vector algebra is applied componentwise. Let X = {Xi, X^) be a 
continuous d-dimensional vector random variable and A = {l,...,d}. We adopt the marginal 



transformation to approximate Laplace margins (Keef et al. , 2012) 



Yi 



log{2FxAXi)} forFx,(X,)<i 
log{2[l-FxAX^)]} i0TFx,{Xi)>\, 



(1) 



where the estimated distribution function Fx^ is obtained from the semi-parametric model of 



Coles and Tawn (1994) 



1 - {1 - FxAux^]{l + ii{x - nxJ/<Ti} , for x > ux,, 



for X < uxi , 



(2) 



where Fxi is the empirical distribution function and uXi is a threshold above which the gen- 
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eralised Pareto distribution with scale parameter and shape parameters o"j and ^j, short-hand 
GP((Ti,^i), (Tj E (0,00), S M, is fitted to the observed values of the excess random variable 
Xi — uXi\Xi > ( [Davison and Smith 1990). The choice of the Laplace transformation is 
motivated by the symmetry of the Laplace distribution that ensures the limiting conditional 



dependence model to be unchanged for negatively dependent variables (Keef et al. 2012). 



2.2 Conditional modelling of extreme values 



The Heffernan and Tawn (2004) conditional dependence model characterises the probabilistic 



behaviour of the conditional vector random variable Y -i\Yi = y, for large y, where Y-i denotes 
the {d— l)-dimensional vector of the transformed variables without the i-th margin. According 



to [Heffernan and Tawn (2004), for each i € A, there exist vector- valued normalising functions, 

> M"'~^, such that for x > 



a 



and 6|j : 



Yi — U > X, : TTT^ < Z 



Yi > ^ exp {—x) G|j (z) as — )• 00, 



(3) 



where the jth marginal distribution of G|j is a non-degenerate distribution function for all 
J G A \ {i} and additionally, the following condition is required such that G|j is well-defined: 

lim Gju (z) = 1 for all j 7^ i 



so there is no mass at +00 but some is allowed at —00, in any margin. Heffernan and Tawn 



( 2004 ) identified the normalising functions are unique up to type, and that for a broad class of 



distributions these functions are all in the parametric family 



aij (x) = aijX and bu{x) = x^^% x>0 



with (Q;|j,/3|j) G [—1,1]"'^"'^ X (— oo,l)''~^. Positive and negative dependence between variables 
Yi,Yj, for i j is given by aj^i > and aj^i < 0, respectively, with aj^i the associated Q|j with 
Yj variable. The strongest form of positive (negative) extremal dependence occurs when aj|j = 1 
(aj|j = —1) and = and is termed as asymptotic positive (negative) dependence, for all 
j 7^ i. Otherwise, variables are termed asymptotically independent. The conditional model of 



Heffernan and Tawn (2004) can be viewed as a multivariate semiparametric regression of Y. 



on Yi, i.e. given Yi > u, for large u we have: 



aux + x^l'Zu lor Yi = X > u (4) 



where Z^^ is a d — 1 dimensional variable with non-zero mean and distribution function 



The original procedure of Heffernan and Tawn (2004) for estimating the vector parameters 



Q;|j and /3|j consists of using pseudo-likelihood methods to jointly estimate the parameters of 
interest. In particular, if we assume that Z^^ has finite vector mean and standard deviations 
<T|j, then the mean and standard deviation of the conditional random variable > u is 

a^^Yi + fj,^i (li)^i' and Yf^^cr^^, respectively. Under the false working assumption that Z|j are 
independent Normal random variables, numerical maximisation of the likelihood of the model 
over the parameter space is required to get parameter estimates (Q;|j, /3|j, <T|j), and G|j is 
estimated nonparametrically by the empirical distribution function of: 

Z\^ = ^ . (5) 

{Y,f\^ 

Once estimates are obtained, standard procedures for inference and extrapolation can be per- 
formed as in Heffernan and Tawn (2004) by implementing Algorithm [Tj As an example, the 



functional P (X € C\Xi > s) can be approximated by repeating steps 1-5, and evaluating the 
estimate as the long run proportion of the generated sample that falls in a set C E M^. As far as 
the confidence intervals of the estimate of any functional are concerned, these are obtained by 
the replication of the three stages of the following bootstrap method: data generation under the 
fitted model, estimation of model parameters and the derivation of an estimate of any derived 
parameters linked to extrapolation. 

Algorithm 1 Sampling Algorithm 
1: Simulate Yi from the Laplace distribution conditional on its exceeding threshold u > 0. 

2: Sample Z|j from G|j independently of Y^. 
3: Obtain Y^i = a\,Yi + {Yif\^Z\i. 

4: Transform Y = iY -i, Yi) to the original scale by using the inverse transformation of ([T]) for 
each margin. 

5: The resulting transformed vector X constitutes a simulated value from the conditional 
distribution of > t^^(n), where denotes the inverse transformation of (nj). 
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2.3 Inference based on Keef et al. (2012) 



Although the efficiency of the model has led to its implementation in a wide range of applications 



including riverflow and rainfall (Keef et al. , 2009), temporal river flow cases (Eastoe and Tawn 



2012), food safety (Paulo et al. , 2006) and finance (Abbas et al. , 2011 ), it was recently discovered 



by Keef et al. ( 2012[ ) that further constraints on the parameter space of the model are required. 
According to Sibuya (1960) and Tiago de Oliveira (1962/63), there are different categorisations 
of extremal dependence between two random variables (Xj, Xj), i.e. asymptotic dependence and 
asymptotic independence measured by the coefficients of tail dependence 



= limP{x, < F-\l-p)\X, > Fr\p)^ 



(6) 



When > {x^j > 0) the variables are termed asymptotically positive (negative) dependent 
and asymptotically independent, otherwise. Taking these measures into consideration, Heffernari] 



and Tawn (2004) omitted the fact that there is stochastic ordering between asymptotically 
independent and dependent models. In particular, let the qth conditional quantile of Yj\Yi = x, 

(2004) model be Vjuiq) = djux + x^^^'^Zju{q), the 



for large x under the 



Heffernan and Tawn 



associated quantile under asymptotic positive dependence y^j(g) = x+z^^-{q) and the associated 
quantile under asymptotic negative dependence yj^^ig) = —x + zT^^ig)- The natural restriction 



yj\i{q)<yj\i{q)<yl\,{q), Vge [0,1], 



(7) 



imposes further constraints on the parameter space of the model which are given by (Theorem 
l.l,|Keef et ar|([20l2l)): 



Case I: either 

aj\i < min |l, 1 - l3jiiZjii{q)v^^\'-\ 1 - v^^^^'^ Zjniq) + v-^z+^.{q)^ . 



7 



or 



l-P,\,Zj\MV'^''^ < ajH < 1 and (l-/3-i){/3.|,z.|,(g)}V(i-/3,|.)(i_a^.|j-/3.K/(i-/3,l.)+^+ (g) > q. 



Case II: either 



0(j\i < min 



{1,1+ l3,\iV^^^^-hj\,{q), 1 + v^^\^-hj\i{q) - v-^zj^.{q)} 



or 



< _«^.|^ < 1 and (l-/3-i)(-/3^.|.z.|.(g))V(i-/3,|.)(i+a.|^)-/3,|./(i-/3,|.)_^- (^^ > q. 



where > u is a value above the maximum observed value of K- so that the constraints are 



imposed only on extrapolations. As far as the selection of q is concerned, Keef et al. (2012) 



found that for both cases the conditions were satisfied for all q if they were each satisfied for 
both = and o = 1. 



3 Estimation of Heffernan and Tawn (2004) model under stochas- 
tic ordering 

3.1 Quantile Ordering Constraints 



In this paper we exploit the same idea for the construction of the parameter space of Heffernan 



and Tawn (2004) model under the assumption of stochastic ordering between two conditional 



random variables. Specifically, let the qth, q £ [0,1], conditional quantile of Yi\Yj = x and 
Yk\Yi = X for large x be yi\j{q) and yk\i{Q), respectively. Under the 



Heffernan and Tawn 



(2004) 



model we have that yi\j{q) = onyx + x^^\^ ziy{q) and yk\i{q) = d^x + Our objective 

is to derive the constraints under which there is stochastic ordering between the conditional 
variables so that the following condition is always satisfied for all x above a level v > u 



yk\i{q) < yi\j{q), Vg g [o, i]. 



(8) 
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The motivation for exploring inequality ([s]) stems from the dose ordering effect in the joint 
region of ALT and TBL. Consider for example the transformed with respect to equation ([T]) 
ALT and TBL and let and y^^((7) be the conditional quantiles of TBL given a large 

level of ALT for dose A and B, respectively. Then under the assumption of liver toxicity, it 
is intuitive to consider the natural ordering of the conditional quantiles < The 

following theorem gives conditions under which two conditional quantiles based on |Heffernan 



and Tawn (2004) satisfy the ordering constraint dSl), for a g G [0, 1]. 



Theorem 1. Let D{x) : [v,oo) — t- M such that D{x) := {aiy — a^i)x + x^'-^^ Ziy{q) — x^''^'zm{q), 
with {aiy,ak\i) e [-1,1?, (A|i> ^ (-oo, l)^ and {zi\j, z^^) e R'^ . For v > u and q € [0,1] 
i) the stationary points (s.p.) of D{x) can be categorised as in Table^ 

Table 2: Conditions for stationary points of D{x). 





number of s.p. 


D'{v) 


s 


D'{s), s G M 







> 


complex/real 


(0,oo) 




1 


< 


complex/real 


(— oo, oo) 




2 


> 


real 


(-00,0) 



where 

and C denotes the set of complex numbers, 

a) The ordering constraint (pi) holds for all x > v if aiy > Qfc|j and either 



1. when D{x) has no s.p.: D{v) > 0, or 

2. when D{x) has one s.p. x^,: inm{D(v), D{x^:)} > 0, or 

3. otherwise: min /^(x*), > 0, where x* and the two s.p. 



Proof 

(i) According to Descarte's rule of signs D'{x) = can have at most two solutions, therefore 
D[x) can have at most two stationary points. Numerical inspection of the function (e.g. for 
ai\j = 0.2, a^i = 0.1,/3;|j = 0.2,/3fc|j = 0.5, Zi\j{q) = 0.6 and zm{q) = 0.6) shows that there can 
be cases where D{x) has two .stationary points. The cases of Table [2] follow from noting that 
D"{s) = is the unique root of D"{x), so that D'{x) has at most one stationary point, i.e. 
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when s > V £ M then s is a s.p. of D'{x), otherwise s is a complex number so that D'{x) is 
monotone for x > v. 



(a) D{x) > for X € [v, oo), implies that lima;_s.oo D{x) can be either or oo, so that: 

ai\j > ak\i. (9) 

Categorising the cases with respect to the number of stationary points of D{x), we have that 
D{x) > 0, for all x > v, if and only if one of the 3 conditions of Theorem^ (ii) holds. □ 

From a statistical perspective, constraints follow from the nature of the D{x) function, i.e. 
one needs to find the stationary points of D{x) so that estimation of Heffernan and Tawn ( |2004 ) 



parameters under the quantile ordering assumption can be carried out. Since the function D'[x) 
is not linear, closed form roots of D'{x) = do not exist. Theorem [T] (i) gives the number of 
s.p. of D[x) together with the intervals that these points can lie, e.g. if D{x) has one s.p. then 
one dimensional root finding is sufficient to estimate the root of D'{x) and if D{x) has two s.p. 
the domain of the function D{x) can be separated into two subintervals (u, s) and (s,oo) both 
in which one dimensional root finding is sufficient to yield estimates of these two s.p.. 

3.2 Inference based on stochastic ordering assumption 



As far as the estimation of the Heffernan and Tawn ( 2004 ) model under stochastic ordering is 



concerned. Theorem [T] provides a set of exclusive cases where each one shows the number of 
stationary points that the function D{x) can have. This provides an automatic way for selecting 
the associated stochastic ordering constraint from Theorem [T] under which maximisation of the 
likelihood of the model is performed. For the required stochastic ordering ([s]) constraint we found 
that the conditions of Theorem [T] were satisfied for all q if they were satisfied for both q = Q 
and 1. To illustrate the feature Figure [T] shows the profile loglikelihood surface for combinations 
of the parameters of the conditional dependence model parameters of TBL given ALT for dose 
A, denoted by a^-^ and under the assumption that the conditional quantile of dose A is 
smaller than the conditional quantile of dose B. Solid lines correspond to the joint g = and 
g = 1 constraints whereas dashed lines correspond to < (7 < 1 constraints. The parameter 
space obtained from the joint g = and q = I constraints is nested in the parameter spaces 
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obtained from < q < 1 constraints. 



- -100 




-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 



Figure 1: Profile log-likelihood surface for dose A parameters (a^^, The solid curves show 
the boundary of the parameter space under the constraints of Theorem [T] when q = and q = 1- 
Dashed curves show the constraints of Theorem 1 when < q < 1, showing these constraints 
are less restrictive than when q = and q = 1- The dot and cross show estimated parameters 
for unconstrained and constrained estimation respectively. 



4 Simulation study of ordering constraints 



The impact of the proposed constraints of Section 3.1 and the Keef et al. (2012) constraints 



is illustrated with a simulation study. We examine the performance of conditional quantile 
estimates using simulated datasets from three bivariate copula models with Laplace marginals. 
Denote by (Yi, Y2) a random variable arising from one of these copula models. We simulate pairs 
of observations conditionally on Yi exceeding a finite threshold u from the exact form of the 
limiting conditional dependence model. Explicitly, we assume that the conditional distribution 
function of I2IY1 > n, for finite n, is equal to the actual limiting distribution function that is 
implied by expression ^ for each of the three copula models of Section 4.1 
that 



I.e., we assume 



1 ~ G211, Yi > u and Z2\i independent of Yi 



(10) 
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with a2|i) (^2\i and G2\i chosen such that expression (Isj) holds. In Section 4.1 we present the 



three copula models along with their Heffernan and Tawn ( 2004 ) limiting representation and 



the simulation design. Algorithm [2] provides the simulation procedure used in this study for the 



bivariate case. In Section 4.2 we discuss the results of the simulation study. 



Algorithm 2 Simulation 
1: Set 1 = 1, (a2|i,/32|i) e (-1,1) X {-oo, 1) and iV e N. 

2: Simulate yij from a Laplace distribution conditional on exceeding a threshold u. 
3: Simulate Z2\ij independently of yij from the true limiting distribution G211 
4: Set y2,i = a2\iyi,i + yij^Z2\ij. 

5: If I < N set 1 = 1 + 1 and go to step 2; otherwise return {yi,y2), where yi 

(yi,i,---,yi,Ar)',?/2 = (y2,i, •••,2/2,Ar)'- 



4.1 Theoretical models 

The three copula models we consider in the simulation study are the bivariate extreme value 
logistic, the inverted logistic and the Gaussian with distribution functions given respectively by 

r{Yi<x,Y2<y)=exp(- [{_ logF(x)}^/'^ + {- log , /^G (0,1], 



P (Yi < x,Y2 < y) =F{X > -x,Y > -y) , where {X, Y) ~ logistic(K), 

F (Fi <x,Y2<y) = $s [^-^ {F{x)] , <^-^ {F{y)]] , 

where F, <I> and <I>s denote the d.f. of the standard Laplace, the standard normal and the 
standard bivariate Gaussian with correlation matrix S = ((Xjj) with cjj^j = 1 and fjjj- = p, 
p G (—1, 1), for i 7^ j. The normalising functions and the residual distribution G2\i of the 
limiting representation ([s]) are summarised in Table ([s]) 

The values of the parameters k, and p used in the simulation study are chosen such that the 
simulated data preserve the stochastic ordering feature. For example, in the logistic copula case, 
dependence increases as the value of k decreases which implies larger joint survival probabilities 
as K decreases. We thus simulate pairs of observations from the exact conditional dependence 



model (10) under k = 0.6 and k = 0.9. For the asymptotically independent models, two pairs 
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Table 3: Heffernan and Tawn (2004) normalising constants a2\i, P2\iy ^''^^ limiting distribution 
G2\i for the bivariate extreme value logistic, inverted logistic and standard bivariate Gaussian 
copula with correlation parameter p > 0. 



Model 


"211 


/32|1 


G2\i(.z) 


logistic 


1 





(1 + exp{-z/«;})''"^ 


inverted logistic 





1 - K 


1 - exp (-Kz^/'^) 


Gaussian 




1/2 


N(0, 2/^2(1 _ ^2)) 



of parameter values are used in the simulation study, i.e., for the inverted logistic copula we 
use K = (0.3,0.7) and k = (1,0.415), and for the Gaussian copula we use p = (0.3,0.7) and 
p = (0,0.5). The second pair of k and p values is used for comparisons between the two 
asymptotically independent models as a key summary measure of extremal tail dependence, the 



coefficient of tail dependence of Ledford and Tawn (1996), between the variables Yi and Y2 is 
the same for these choices. 



The conditional quantile estimates are obtained from the original Heffernan and Tawn (2004) 



model, the constrained model of Keef et al. (2012) and the constrained model described in 



Section 3.1 We refer to these models as HT, AD and SO, respectively. On this note, the 
stochastic ordering constraints are imposed on pairs of observations with different dependence 
parameters. The performance of the estimates is assessed with the Monte Carlo estimate of the 
root mean square error. To be specific let y{q) and y{q) be the true conditional quantile and its 
model-based estimate. The Monte Carlo estimate of the root mean square error is 



RMSE(g) 



^ m 

— y]{yi(9) -y(9)}^ 



where yi{q), ■■■,ymiQ) denotes a Monte Carlo sample of the conditional quantile estimates 
and each estimate is obtained from a simulated sample of G N pairs of observations, 
{{yi,y2) '■ Ui.i > u,i = l,...,N}. For all three models, the root mean square error is estimated 
over a grid values of g G [0, 1] and three conditioning levels, i.e., xo.95, xo.99 and xo.999, where Xp 
denotes the p-th quantile of the standard Laplace distribution. Comparisons are made on the 
basis of ratios of RMSEs between estimates from different models. For the constrained models 
we tabulate the percentage of the Monte Carlo samples where estimates changed with respect 
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to the original Heffernan and Tawn (2004) model. 



4.2 Results of simulation study 

Table |4] shows the percentage of the estimates that changed with respect to one of the three 
reference models (HT, AD, SO). The imposition of constraints to the parameter space of the 
HT model alters estimates particularly in the asymptotically independent models and less in 
the asymptotically dependent model. In particular, the larger changes occur when variables are 
highly dependent except for the logistic model. Regarding the logistic model, the percentage 
of changes in the first two rows is relatively small compared to the other models. This feature 
is explained by the greater concentration of the maximum likelihood estimates around the true 



value, which lies in the feasible set of the constrained space (Keef et al. , 2012 ). Additionally, small 
changes occur within the asymptotically independent models especially when the variables do 
not possess strong dependence. For the second pair of parameter values for the inverted logistic 
and Gaussian copulas, the changes in parameter estimates do not occur at a similar rate when 
dependence between variables is present {k = 0.415, p = 0.5). We therefore conclude that the 
constraints induced by the AD and SO models are not only related to the level of dependence 
but to the dependence structure as well. Table [5] shows the ratio of the Monte Carlo root mean 



Table 4: Percentage of estimates that changed with respect to a reference model. First row: 
percentage of AD estimates different from the HT estimates, second row: percentage of SO 
estimates different from the HT estimates, third row: percentage of SO estimates different from 



AD estimates. Columns show the corresponding copula model from Section 4.1 used in the 
simulation. 





logistic copula 


inverted logistic copula 


Gaussian copula 


K or p 


0.6 


0.9 


0.3 0.7 


0.415 1 


0.7 0.3 


0.5 


AD-HT 


29% 


27% 


63% 10% 


41% 0.3% 


36% 0% 


6% 1% 


SO-HT 


54% 


56% 


77% 45% 


68% 47% 


42% 10% 


30% 25% 


SO- AD 


33% 


46% 


43% 44% 


47% 47% 


9% 9% 


24% 25% 



square error, of the conditional quantile estimates obtained from the three copula models. An 
increase in efficiency under the imposition of the constrained models AD and SO is observed 
for nearly all conditional quantile estimates in the asymptotically independent models. The 
highest reduction in RMSE is achieved by the SO model in the inverted logistic copula, a 
feature which is also consistent with the higher percentage of change in estimates as shown 
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in Table |4j The conclusion for the asymptotically independent models is that the efficiency 
of the conditional quantile estimates is, in decreasing order, SO, AD and HT. Regarding the 
asymptotically dependent logistic copula, constrained models appear to be less efficient than 
the HT model and the efficiency of the conditional quantile estimates is, in decreasing order, 
HT, AD and SO. 

Table 5: Ratio of the Monte Carlo root mean square error of the conditional quantile estimates 
y2|i(9) obtained from the HT, AD and SO models. Results are reported for q = 0.2, 0.5, 0.8 and 
three conditioning levels xo.95, 2:0.99 and xo.999. The value Xp here denotes the p-th quantile of 
the standard Laplace distribution. 



I logistic copula | inverted logistic copula | Gaussian copula 

Korp I 06 09 r03 07 0.415 1 I 0.7 03 05 



q = 0.2 



AD/HT 


xo.95 1-10 1.02 
X0.99 1.00 1.08 
X0.999 1-10 1.04 


0.99 1.00 0.98 1.00 
0.99 0.97 0.98 0.98 
1.00 0.96 0.97 0.98 


1.00 1.00 0.99 0.99 
1.00 1.00 1.00 1.00 
0.99 1.00 0.99 0.99 


SO/HT 


xo.95 1.20 0.97 
X0.99 0.99 1.24 
X0.999 0.96 1.12 


0.93 0.99 0.97 0.96 
0.95 0.97 0.94 0.96 
0.96 0.94 0.94 0.98 


0.98 0.99 0.99 1.00 
0.97 0.99 0.99 1.00 
0.97 0.99 0.97 0.97 


SO/ AD 


xo.95 1-10 0.94 
X0.99 0.96 1.14 
.X0.999 1.01 1.07 


0.93 0.99 0.99 0.96 
0.95 0.99 0.95 0.97 
0.96 0.98 0.96 0.99 


0.98 0.99 0.99 1.00 
0.97 0.99 0.99 1.00 
0.98 0.99 0.97 0.97 


q = 0.5 


AD/HT 


Xo.95 0.99 1.05 
X0.99 1.09 1.08 
X0.999 1-11 1-08 


1.02 1.02 0.98 1.00 
1.01 0.96 0.98 0.98 
1.00 0.89 0.94 0.95 


1.00 1.00 0.99 1.00 

1.01 1.00 1.00 1.00 
1.00 1.00 0.99 0.97 


SO/HT 


xo.95 0.78 1.43 
X0.99 1.09 1.29 
X0.999 1.15 1.27 


0.85 1.02 0.89 0.98 
0.79 0.96 0.82 0.93 
0.84 0.85 0.82 0.92 


0.94 0.99 0.94 1.01 
0.90 0.98 0.93 1.00 
0.86 0.98 0.90 0.90 


SO/ AD 


xo.95 0.78 1.35 
X0.99 0.99 1.19 
X0.999 1.03 1.17 


0.84 0.99 0.89 0.97 
0.78 0.99 0.83 0.94 
0.84 0.95 0.86 0.96 


0.94 0.99 0.94 1.01 
0.89 0.98 0.93 1.00 
0.86 0.98 0.91 092 


q = 0.8 


AD/HT 


xo.95 1.04 1.05 
X0.99 1.10 1.08 
X0.999 1-14 1.09 


1.02 1.03 1.00 0.99 

1.03 0.98 0.99 0.99 
0.81 0.86 0.94 0.93 


1.00 1.00 0.99 1.00 

1.01 1.00 1.00 0.99 
1.00 1.00 0.99 0.95 


SO/HT 


xo.95 0.97 1.25 
X0.99 1.12 1.29 
X0.999 1.18 1.31 


0.84 1.04 0.84 0.99 
0.73 0.99 0.77 0.95 
0.71 0.82 0.75 0.88 


0.91 0.99 0.92 1.02 
0.87 0.98 0.90 0.99 
0.81 0.97 0.88 0.88 


SO/ AD 


xo.95 0.93 1.19 
X0.99 1-01 1.19 
X0.999 1-04 1.19 


0.81 1.00 0.84 1.00 
0.71 1.00 0.78 0.95 
0.74 0.94 0.79 0.94 


0.91 0.99 0.92 1.01 
0.86 0.98 0.90 1.00 
0.80 0.97 0.88 0.92 
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5 Application: drug-induced liver injury 
5.1 Preprocessing and outline of analysis 

The data that we consider in this study relates to a sample of 606 patients that were issued a 
drug linked to liver injury in a randomised, parallel group, double blind Phase 3 clinical study. 
ALT and TBL measurements were collected from all patients at baseline (prior to treatment) 
and post-baseline (after 6 weeks of treatment) periods. Let V^^ and V-^p be the i-th. baseline 
and post-baseline laboratory variable respectively, measured at dose j = A, B,C and D. We 
use i = 1, 2 to denote the ALT and TBL, respectively. 



Instead of working with the raw data, the Box and Cox (1964) transformation is applied ini- 



tially to stabilise the heterogeneity observed in the samples. For this dataset we apply the log- 
transformation and we denote the transformed data by W- ^ = log(y/^) and W- p = log{V/p). 
Consequently, we use a robust linear regression model of the log-post-baseline on the log-baseline 
variable to adjust for the baseline effect, i.e., in its simplest form, the robust linear regression 
of W- p on W^p admits 



Wlp = li + SiWlp + Xl, i = l,2, j 



A,...,D, 



(11) 



where (7^, 5^) E and Xj is a zero mean error random variable. Here we use median quantile 



regression (Koenker and Bassett, 1978) which is equivalent to assuming that the error random 



variable Xj follows the Laplace distribution with zero location constant scale parameters rtYu 



and Moyeed 2001 ). Our approach is based on the basic model structure of Southworth and Hef- 



fernan 



(20126), i.e., the extremal dependence of {Xl,X2) is estimated from the 



Heffernan and 



Tawn 



(2004) conditional dependence model whereas the log-baseline variables p and W2 p 



are modelled independently for each dose. Under the assumption of independence between X^ 
and W- p simulated samples of the post-baseline variables can be generated. The exact proce- 
dure of the simulation is straightforward, i.e., residual and baseline samples are generated from 
their models and are combined in equation ( |ll| , with and dj replaced by their corresponding 
maximum likelihood estimates, to produce simulated samples for the log-post-baseline variable 
W-p. The simulated sample is then back-transformed to its original scale using the inverse 
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Box-Cox transformation. 



The key differences between our modelling procedure and Southworth and Heffernan (20126) are 



related to the modelling of the baseline and the estimation of the conditional dependence model 
parameters. Firstly, for each baseline variable we implement the univariate semi-parametric 
model of Coles and Tawn (1994) as described in Section 2.1 by equation ^ whereas Southworth 



and Heffernan (20126) use the empirical distribution function. Our motivation for modelling the 



tail of the baseline variable stems from the fact that it is likely to observe higher baseline ALT 
and TBL in the population (post-marketing period) than in the clinical trial (pre-marketing pe- 
riod). Therefore, tail modelling of the baseline is key to the simulation process as it incorporates 
a natural source of extremity through model-based extrapolation. Results from the univariate 



analysis are not presented in this paper but similar analyses can be found in Southworth and 



Heffernan (2012 a) and Papastathopoulos and Tawn (2012). As far as the estimation of the 



conditional dependence model's parameters is concerned, we impose the stochastic ordering 
constraints developed in Section [3] and compare the results with the unconstrained estimates 
obtained from the HT model. 



In Section 5.2 we illustrate the effect of the stochastic ordering constraints via estimates of 
conditional quantiles for all doses and in Section [5. 3| we proceed to the prediction of the prob- 
ability of DILI by simulating post-baseline laboratory data of hypothetical populations of size 
200000 using the fitted marginal and conditional dependence models. The assessment of the 
uncertainty of the estimates of extreme quantities is performed via the bootstrap procedure. 



5.2 Dependence modelling 

Let and Y,^ be the transformed, with respect to equation ([l|, residuals and for each 
dose j. Figure [2] shows the bivariate scatterplots of against Y-^ for all dose levels. The 
dependence between the residual ALT and TBL variables appears to be very weak for all dose 
levels and a direct conclusion regarding the stochastic ordering effect cannot be made on the 
basis of Figure [2j A lack of ordering appears from the estimated conditional quantiles of TBL 
given ALT from the HT model as shown in Figure|3]in the standard Laplace scale. The estimates 
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of the median conditional quantiles from the HT model are ordered above approximately the 
conditioning level 4 whereas the minimum and maximum conditional quantile estimates exhibit 
a lack of ordering for the majority of the conditioning levels. The effect of constraining the 




I- 



2 4 

rD 



Yt 



Figure 2: Scatterplots of Y2 against Y"/, for j = A, D 

parameter space to impose the stochastic ordering assumption between all dose levels is shown 
in Figure [S] via the conditional quantile estimates obtained from the SO model. For the SO 
model we selected v to be 5, the 99.7% quantile of the Laplace distribution. The imposition of 
the ordering constraints induces changes in all conditional quantile estimates which satisfy the 
ordering assumption above the conditioning level 5. The most important change in the quantile 
estimates is observed for dose A which are considerably smaller than the HT estimates, when 
q = l. 



5.3 Prediction 

In this section, the main focus is placed on the prediction of joint elevations of ALT and TBL. As 



stated by FDA (2008) and mentioned earlier in the Introduction, DILI is associated with ALT 



and TBL exceeding the SxULNalt and 2xULNtbl respectively. For ALT, the ULN is taken 
to be 36 units/litre and for TBL is 21 ^umol/litre. Let be the joint survival probability 
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q = q = 0.5 9 = 1 




2468 10 2468 10 2468 10 



q = q = 0.5 Q = ^ 




2468 10 2468 10 2468 10 



Figure 3: Conditional quantile estimates y2\i{q) of Y2\Y-^ = x under HT (first row) and SO 

models (second row). The line types correspond to - • - • - dose A, dose 

dose C, dose D. First, second and third column show the minimum [q = 0), median 

{q = 0.5) and maximum [q = 1) conditional quantiles, respectively. 
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of {ALT > X n TBL > y} at dose level j = A,B,C or D, i.e., 



p>ix,y) ■.= F{VIp> x,V^p> y), x,y G 



(12) 



To estimate the survival probability (12) we follow the approach of Southworth and Heffernan 



(20126), also mentioned earlier in Section 5.1 and simulate N = 200000 post-baseline samples. 



For each dose level, baseline samples are generated from the semi-parametric model ^ and 



are subsequently combined with N generated residual samples from the SO constrained Heffer 



nan and Tawn 



(2004) model in equation (11), with and 6^ replaced by their corresponding 



maximum likelihood estimates, to produce simulated samples for the log-post-baseline variable 
W- p. The simulated sample is then back-transformed to its original scale and the survival prob- 



ability (12) is estimated empirically. To assess the uncertainty of the estimates, this procedure 



is repeated 1000 times and 95% equal-tail confidence intervals are obtained from the bootstrap 
distribution of each estimate. 



Figure [4] shows the estimated survival probabilities p'{x,y) for x = 3 x ULNalt and variable 
y. For comparisons, estimates are reported from the SO and HT models. The imposition of 
the constraints induces changes in all estimates. In particular, the survival probability esti- 
mates from the SO model are lower than the HT model for all doses, especially in the region 
{y : y < 30}. This behaviour also implies changes in the upper tail and in the joint region of 
DILI, i.e., when x = 3 x ULNalt and x = 2 x ULNtbl- 



To assess the effect of the stochastic ordering constraints in the joint upper tail of ALT and 
TBL we study in Figure [s] the Monte Carlo estimates of extreme quantiles y^ defined by 

2/^ : F{VIp > 3 X ULNalt, V^p > y^) = l/T, 

for large T. The level y^ is more commonly referred to as the return level associated with 
patient return period T and is expected to be exceeded, jointly with {U/p > 3 x ULNalt}, 
once in T patients. 
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Figure 4: Estimated survival probabilities p'{x,y) for all doses with x = 3 x ULNalt and 
variable y. Solid and dashed lines correspond the point-wise estimates and their 95% confidence 
intervals, respectively. Black and grey colour shows estimates from the SO model and HT 
models, respectively. 




7700 44000 260000 1200000 1.8e+07 7700 44000 260000 1200000 1.8e+07 

Patient return period Patient return period 



Figure 5: Monte Carlo return level estimates of return levels plotted against the patient 
return period T on the logarithmic scale. Estimates are obtained from the SO (left) and HT 
(right) models and the line types are as in Figure ^ and correspond to the dose levels. 
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The figure indicates a lack of ordering in the estimates obtained from the HT model which 
stems from the behaviour of the return level estimates of dose B that are unexpectedly higher 
than doses C and D for large patient return periods. On the other hand, the return level esti- 
mates obtained from the SO model satisfy the ordering constraints and are larger than the HT 
estimates. On this note, we also found smaller return level estimates than those predicted by 
the HT and SO models when samples of the baseline variables where generated via resampling 
from the observed samples and not from the semi-parametric model ([2]). This is attributed to 
the fact that the resampling approach does not allow for extrapolation in the baseline variables. 

We view this lack of ordering as a by-product of sampling variability that arises in small samples 
and thus the SO model has corrected for this. 
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